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ABSTRACT 

A  continuous  identification  of  parameters  is  performed 
on  a  simulated  fast  breeder  nuclear  reactor  system  using 
hybrid  computation  and  applying  techniques  of  statistical 
regression  analysis  and  exponentially-mapped-past  functions. 
Output  states  which  are  not  directly  measurable  are  estimated 
by  use  of  a  Kalman  filter. 

The  method  developed  in  this  study  is  applied  to  a  nu- 
merical example  which  demonstrates  that  unknown  parameters 
can  be  identified  within  3%  of  their  actual  value,  with  sig- 
nal  noise  ratios  as  low  as  10:1  in  the  measured  states.   The 
example  also  demonstrates  that  convergence  occurs  in  a  rea- 
sonably short  time. 

An  executive  software  routine  was  written  in  FORTRAN  IV 
for  the  XDS-9300  digital  computer  and  was  applied  both  for 
the  control  of  the  process  and  for  establishing  the  opera- 
tional performance  and  accuracy  of  the  Comcor  Ci-5000  analog 
computer  that  simulated  the  nuclear  reactor  system. 
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I.   INTRODUCTION 

Fast  breeder  nuclear  reactors  such  as  the  recently 
proposed  design  described  in  [1]  do  not  have  inherent  load 
following  ability  and  therefore  require  an  automatic  con- 
troller to  handle  load  changes.   It  is  important  in  design- 
ing for  stable  operation  of  the  automatic  control  system  to 
know  the  value  of  the  parameters  of  the  plant  within  a  rea- 
sonable degree  of  accuracy.   These  parameters  change  due  to 
such  factors  as  temperature,  fuel  depletion,  fuel  loading 
and  power  level,  [2],  making  it  desirable  to  have  continuous 
identification  of  parameters. 

The  continuous  identification  of  parameters  of  dynamic 
systems  has  received  intensive  attention  during  the  last 
few  years  [3],  [4].   References  [5]  and  [6]  among  others, 
have  shown  successful  results  by  applying  gradient  techniques 
and  quasilinearization,  respectively.   These  methods  are  cum- 
bersome for  nonlinear  systems  and  in  particular  for  dealing 
with  a  nonlinear  space-dependent  reactor  system,  [6],  where 
the  requirement  of  continuous  identification  cannot  be  met 
since  the  in-core  sensors  which  are  necessary  for  certain 
state  measurements  have  not  yet  been  developed1 . 

The  implementation  for  continuous  parameter  identifica- 
tion is  difficult  and  complex.   For  example,  many  gradient 


i 

The  procedure  available  consists  of  the  exposure  of 
a  matallic  foil  by  immersion  and  takes  enough  time  (between 
one  and  two  hours)  to  eliminate  any  possibility  for  con- 
tinuous identification  of  parameters. 


techniques,  even  when  they  represent  in  some  cases  very  ele- 
gant ways  of  approaching  the  hill-climbing  search  procedure 
are  limited  by  practical  considerations  in  their  application. 
To  illustrate  this,  reference  is  made  to  a  particular  gradi- 
ent (steepest  descent)  technique  which  uses  an  error  crite- 
rion and  sensitivity  equations,  and  which,  among  other 
requirements,  has  the  need  for  a  forcing  function.   A  block 
diagram  of  the  algorithm  is  shown  in  Figure  1.   In  this 
method  it  is  required  to  derive,  implement  and  solve  the 
sensitivity  equations  in  order  to  obtain  the  value  of  the 
gradient  for  the  specified  cost  function.   In  order  to  iden- 
tify n  parameters  the  number  of  analog  components ,  using 
standard  methods,  must  be  multiplied  by  at  least  n  *  ,  since 
for  a  linear  case  these  equations  take  the  same  form  of  the 
model,  [7].   In  addition,  for  nonlinear  systems  the  sensi- 
tivity equations  have  more  complicated  forms  which  do  not 
permit  the  use  of  techniques  such  as  multiplexing,  [7].   To 
apply  these  methods  the  system  must  be  initialized  with  some 
previous  information,  usually  an  educated  guess  of  the  pa- 
rameter values,  and  then  the  process  must  be  iteratively 
reinitialized  until  some  criterion  is  satisfied.  Since  the 
parameter  correction  is  a  function  of  the  gradient  and 
since  this  gradient  approaches  zero  as  the  parameter  value 
approaches  the  actual  value  (solution) ,  the  procedure  will 
converge  in  a  reasonable  time  only  if  a  means  is  provided 


More  recent  techniques,  [8],  allow  for  simplification, 
when  the  model  can  be  expressed  in  canonical  form. 
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to  speed  up  the  convergence  in  the  vicinity  of  the  minimum 
of  the  cost  function. 

Other  techniques,  [9]  and  [10],  which  are  used  in  non- 
linear systems  suffer  from  more  or  less  the  same  problems. 
These  difficulties  indicate  the  need  for  an  identification 
system  that  will  operation  without  supervision,  require 
neither  educated  guesses  nor  initialization,  be  very  reli- 
able, and  use  simple  implementation  that  avoids  the  trouble- 
some iterative  initializing  process. 

The  procedure  developed  in  this  thesis  makes  use  of 
the  method  of  statistical  regression  analysis,  applied  pre- 
viously with  success  to  linear  systems,  [4],  and  confined 
for  a  long  time  to  the  analysis  of  statistical  populations, 
[11]  and  [12].   The  method  of  steepest  descent  is  applied 
to  the  classical  regression  analysis  in  order  to  speed  up 
the  convergence.   To  avoid  the  need  for  reinitialization, 
the  properties  of  exponentially-mapped=past  statistical 
variables  is  used,  [13]. 

Although  the  procedure  is  applied  here  to  a  simplified 
model  of  a  fast  breeder  nuclear  reactor  system,  the  method 
is  general  enough  to  be  applicable  to  many  other  automatic 
control  systems  that  are  represented  by  more  complex  models 

For  the  particular  nuclear  reactor  system  used  as  the 
example  application,  the  variable  representing  the  delayed 
neutron  concentration  is  a  nonmeasurable  state  and  in  order 
to  have  its  value  for  the  identification  system  a  Kalman 
filter,  [14]  and  [15],  was  used. 
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The  unknown  parameters  were  identified  within  3%  of 
their  actual  values.   The  example  demonstrates  that  con- 
vergence occurs  and  that  the  time  of  convergence  compares 
favorably  with  current  methods,  [16].   The  simulation  re- 
sults shown  in  this  thesis  were  obtained  using  zero  as  the 
initial  values  for  the  parameters;  however,  tests  have 
shown  that  the  convergence  time  can  be  greatly  improved  if 
an  educated  guess  is  accepted. 

An  executive  software  routine  was  written  in  FORTRAN 
IV  for  the  XDS-9300  digital  computer  and  was  applied  both 
for  the  control  of  the  process  and  for  establishing  the 
operational  performance  and  accuracy  of  the  Comcor-5000 
analog  computer  that  simulated  the  nuclear  reactor  system 
and  was  used  to  implement  the  analog  part  of  the  identifi- 
cation procedure. 
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II.   REACTOR  DYNAMICS 


The  point  model  of  the  nuclear  reactor  kinetic  equa- 

2  3  5 

tions  for  U    fuel  can  be  found  in  several  references, 

[16]-[19]. 

These  equations  are 

6 

dn   _    6K      n        0      „         V      ,    r  M. 

ttt  =   irx     n-inrn+Z_l     A.C.  (1) 


and 


It  "  V?    "  I*  i^i 


dC    3- 


T    =  7T3T  n  "  XiCi  i  =  1,  .  .  .  ,6.  (2) 


dt   I*   ll    l  i 

The  left-hand  term  in  equation  (1)  represents  the  rate 
of  change  of  the  neutron  density;  the  first  term  on  the 
right  side  of  equation  (1)  represents  the  rate  of  production 
of  prompt  neutrons  ,  the  second  term  the  rate  of  absortion 
plus  leakage  and  the  last  one,  the  production  rate  of  de- 
layed neutrons.   The  variable  6K,  the  reactivity  input,  can 

be  expressed  as 

I 

6K  =  p   -  L,       a.T. . 
e    j-1    3  ^ 

Applying  this  substitution,  equations  (1)  and  (2)  become 

I  6 


dn 


lPe 


I  Vjp"  ^n  +  £  "i'i  (3) 


dC.    3- 

-±   =  yi  n  -  X.C.  (4) 

dt    l*  l  l  v 
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where   T    is    defined  by 


dT.  / 


-    1 
o  / 


M'jkh  -Tk]    J-!.-.*- 


Since  u.  and  3-  are  fractional  quantities  it  follows  from 
j       i  n 

their  definitions  that 


and 


,t. 


=  B 


The  initial  conditions  for  equations  (3)  thru  (5)  are 


n(0)  =  n 

v  J  o 

B-  n 

c-  (0)  =  c.   =  nl     ° 
i v  ^     io    £*  A  . 


i  =  1,  ...  ,6 


Tj(0)  =  0  j  =  1,...  ,£. 

For  the  purpose  of  this  study  a  single- time-constant 
feedback  reactor  with  a  single  group  of  delayed  neutron 
precursors  will  be  used;  the  following  equations  define  the 
system 


JY  =   Pe  "  aTl  <r£    or  n  +  XC 


2,1 


Jl! 


(6) 


dC 
dt 


jt  =  tttt  n  -  XC 


dT 

dt 


£* 


=  y 


n 


n 


aT 


(7) 
(8) 
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The  initial  conditions  are 

pe(0)  =  0 
n(0)  =  n 

cw  ■  co    mr  "  . 

T(0)  =  0. 

By  changing  variables  as  indicated  in  [16]  ,  (see  Ap 
pendix  A)  a  more  suitable  representation  is 


N  -  A   [(K   -  aT  -  1)  N  +  D]  (9) 


D  =  X  (N  -  D)  (10) 

T  =  Y  (N  -  1  -  bT),  (11) 


where 


xt    rip.    C        v  pe    3    c  1n3       ot       y 

N  -  — ■-.  D  =  jr-   ,  Ke  =  -j   ,  jg-  =  6x10   ,  a  =  ^  ,  Y  =  H. 


b  =  -f  , 

y 

and  the  initial  conditions  are 

N(0)  =  1  T(0)  =  0 

D(0)  =  1  Ke(0)  =  0. 

Experience  with  present  day  design  of  nuclear  reactors 
indicates  that  these  systems  are  stable  under  normal  opera- 
ting conditions.   The  large  breeder  reactors  are  still  in 
the  design  stage  and  a  computer  simulation  analysis  has 
shown  that  these  plants  will  be  stable,  [1] .   Nonetheless, 
a  stability  analysis  applying  the  method  of  Liapunov  has 
been  performed  for  this  model  and  is  included  in  Appendix  B 

16 


The  model  defined  by  equations  (9)  thru  (11)  was  im- 
plemented, but  acceptable  accuracy  could  not  be  obtained. 
The  large  value  of  8/£*  required  a  large  scaling  factor  and 
also  a  large  gain  for  the  state  variable  N.   Because  of 
these  t\co  reasons  the  voltage  that  represented  the  state 
variable  N  become  too  small  and  too  close  to  the  uncertainty 
level  (fourth  digit)  of  the  analog  computer. 

From  simulation  it  is  observed  that  N  is  approximately 
equal  to  zero.   This  and  the  fact  that  B/&*  is  very  large 
implies  that  equation  (12)  is  a  satisfactory  approximation 
for  equation  (9).   This  further  simplification  will,  of 
course,  introduce  some  inaccuracies  since  it  is  assumed,  in 
fact,  that  N  is  zero;  however,  it  has  been  shown,  [16],  that 
the  largest  error  introduced  with  this  assumption  is  less 
than  2%. 

With  this  simplification  the  model  is 

N  =  -  D/(Ke  -  aT  -  1)  (12) 

D  =    (N  -  D)  (13) 

T  =  T(n  -  1)  -  byT.  (14) 

In  state  variable  form  where 

Xj  =  D   ,   x2  =  T 

x1    =    X  (x3  -  xj  (15) 

x2  =  y  (x3  -  1)  -  byx2  (16) 

and  the  variable  x3  is  defined  as 

x,  =  -  x,  /(K   -  ax2  -  1).  (17) 
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For  a  numerical  example,  the  known  constants  are  as- 
sumed to  be  A.  =  0.4  and  y   =  1.0  and  the  parameters  to  be 
identified  are  "a"  and  "b".   Using  these  values  the  equa- 
tions are 

i,  =  0.4  (x3  -  xj  (18) 

x2  =  x3  -  1  -  bx2  (19) 

and  the  variable  x3  is  defined  in  equation  (17) .   The  analog 
computer  equations  using  the  scaled  variables 

Xj  =  40xx 
x2  =  40x2 


and 


are  then 


x  3    Zx  3 


x,  =  8x   -  0.4x  (20) 

•  3  1 

x2  =  20x3  -  40  -  bx2  (21) 

where  the  variable  x   is  defined  by 


2*! 
X3    "  40K   -  ax,  -  40  ^22) 


Figure  2  shows  the  implementation  of  the  above  equa- 
tions and  the  settings  of  the  potentionmeters  for  this 
scaling  are  shown  in  Table   I. 
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TABLE  I 


POTENTIOMETER  SETTINGS 


Pot  Value 

P001  0.200 

P002  0.800 

P003  0.500 

P004  0.500 

P005  0.200 

P006  0.400 

P010  0.200 

P011  0.600 


Pot 

Value 

P012 

0.200 

P015 

0.600 

P017 

0.600 

P021 

0.600 

P051 

0.600 

P055 

0.600 

P057 

0.600 

20 


III.   REGRESSION  ANALYSIS 

Several  attempts  were  made  in  this  study  to  achieve 
the  desired  goals  for  the  identification  system  as  stated 
in  Section  I.   Arrising  in  all  of  these  attempts  was  the 
problem  stated  by  R.  C.  K.  Lee,  [14],  and  illustrated  in 
Figure  3,  which  is  related  to  the  closed-loop  formed  by 
the  estimator  and  parameter  identification  system.   This 
problem  is  difficult  because  the  overall  system  is  nonlin- 
ear even  for  linear  plants.   Because  convergence  of  this 
loop  was  not  proved,  and  because  estimation  and  identifica- 
tion of  parameters  for  nonlinear  cases  depends  strongly  on 
the  initial  value2 ,  it  was  decided  to  approach  the  problem 
in  the  following  alternative  way. 

The  system  is  nonlinear  as  defined  by  equations  (9) 
thru  (11).   However,  the  simplification  performed  on  the 
model  that  leads  to  equations  (18)  and  (19)  yields  a  second- 
order  linear  model  with  the  variable  x3  representing  an  in- 
put.  Since  it  is  possible  to  measure  N  and  T,  which  are 
represented  in  the  simplified  model  by  the  variable  x3  and 
the  state  variable  x2  ,  respectively,  then  the  non-measurable 
state  D,  represented  by  x  ,  can  be  easily  estimated.   In 
order  to  make  the  simulation  realistic,  the  variable  x3  is 
treated  as  a  non-deterministic  forcing  function.   This  func- 
tion is  generated  by  applying  equation  (22)  and  the  resulting 


2 

It  is  possible  for  this  kind  of  system  to  reach  a 
steady-state  value  different  from  the  actual  value  (this  is 
known  as  a  state  of  nonlinear  stability) . 
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value  is  corrupted  with  some  Gaussian  noise  of  mean  zero 
obtained  from  the  low-frequency  Gaussian  noise  generator 
in  the  Comcor-5000  analog  computer. 

This  simplification  has  the  advantage  that  it  opens 
the  estimator-identification  loop,  Figure  4,  because  the 
only  estimated  state  (xx)  is  neither  a  function  of  the  un- 
known parameter  "a",  nor  a  function  of  the  unknown  parameter 
"b". 

The  technique  used  in  this  study,  known  as  statistical 
regression  analysis,  consists  basically  of  finding  the  best 
fit  to  a  set  of  data  using  a  least- squares  criterion.   It 
has  been  shown,  [4]  and  [10],  how  the  steepest  descent 
method  can  be  used  on  the  analog  computer  for  parameter  op- 
timization and  model  building,  and  therefore  this  method 
was  applied  in  the  derivation  of  the  regression  equations 
(27)  and  (30) . 

In  order  to  perform  the  regression  analysis  it  is  re- 
quired to  define  S,  a  cost  function.   The  cost  function, 
for  simplicity,  is  usually  chosen  to  be  the  integral  of 
the  square  of  the  error,  which  is  defined  by  a  relationship 
that  contains  the  parameter  to  be  identified.   The  second 
step  is  to  find  the  rate  of  change  of  S  with  respect  to  a 
new  independent  variable3,  t,  and  to  minimize  the  cost  func- 
tion using  the  method  of  steepest  descent.   The  details  of 
these  methods  are  developed  in  the  following  sections. 


3 

The  variable  t  represents  computer  time  and  it  changes 
more  rapidly  than  real  time. 
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A.    IDENTIFICATION  PROCEDURE  FOR  THE  UNKNOWN  PARAMETER  B 

.From  equation  (14)  an  error  function  z1    can  be  defined 
as 

Ej  =  T  -  y(N-l-bT)  .         (23) 

and   a   cost   function   S2    as 


t  t 

I  [T-T(N-l-bT)]2 

o 


z  z 

i  ^  J   e?  dt  =  Jli 


b  is  an  estimate  of  the  parameter  b. 

The  partial  derivative  of  Si  with  respect  to  time  is 

3S1    dS1         a£ 

9t     3b    3x 

Using  the  method  of  steepest  descent  applied  directly 
to  the  objective  function  it  is  required,  [4],  to  make 

3b    8Si 


3t     3b 

This  ensures  that  as  x  increases  the  estimate  b  changes 
in  such  a  way  that  Si  decreases.   When  the  estimate  b  makes 
3b/3i  =  0  the  minimum  of  the  cost  function  has  been  reached. 

Taking  the  partial  derivative  of  equation  (24)  with 
respect  to  b  gives 


J  li~ 


3Sj 

-x-  =  2     [T- y(N-l-bT)]y  T  dt 
3b 


and,  therefore 
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r  t 


3b 

9x 


=  -2Y 


fTtdt-YjNTdt  +  y  f  Tdt+b   f 

-  0  0  '0  0 


T  dt 


(25) 


Letting 


X,  =  D 

x2  =  T 
x3  =  N 


(26) 


and  substituting  (26)  into  equation  (25)  gives  the  final  ex- 
pression 


9b_ 

9x 


■    2Y 


t  L  t  t 

Y     I      x2x2dt-Y  I   x2x3dt  +  Y/   x2dt+b     I  x2di 

-  „  0  0  0  0 


(27) 


Since  x2  and  x3  are  both  measurable  and  continuously 
available,  the  identification  of  the  unknown  parameter  "b" 
could  also  be  done  continuously  and  the  input  to  the  iden- 
tification system  does  not  require  any  digital  processing. 
Due  to  a  lack  of  analog  multipliers,  however,  all  multipli- 
cations were  performed  in  the  digital  computer.   Figure  5 
shows  the  implementation  diagram  for  equation  (27) .   In 
Section  III-D  a  modification  to  the  procedure  is  introduced 
in  order  to  avoid  the  need  for  reinitialization;  this  modi- 
fication uses  the  properties  of  exponentially-mapped-past 
averages. 

B.    IDENTIFICATION  PROCEDURE  FOR  THE  UNKNOWN  PARAMETER  A 

The  parameter  "a"  can  be  identified  using  equation  (12) 
Since  D  is  a  non-measurable  state  and  it  is  needed  in  the 
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regression  equation,  a  Kalman  filter,  as  described  in  Sec- 
tion JII-C,  was  used  to  obtain  the  estimate  of  the  variable 
D. 

Taking  equation  (12) 


N  =  -  D/(Ke-aT-l) 
and  defining  an  error  criterion  e   as 


z2    £  D  +  N(Ke-aT-l) , 


leads  to  the  cost  function 

t 
S. 


(28) 


/ 


e2  dt. 

2 


(29) 


a  is  an  estimate  of  the  parameter  "a".  The  steepest  descent 
algorithm  can  be  applied  to  obtain  the  differential  equation 
that  defines  the  regression  analysis  for  the  unknown  parame- 
ter "a". 

Using  (26)  gives 


and 


£2  =  txi  +  x3  (Ke"a  x2  "  !)]2  » 


dS2  d 


da        3a 


J    [x,  +  x3  (Ke-a  x2  -  l)]2  dt 


=  2j  [x1+x3  (Ke-2x2-1)](-x2x3)dt 


Therefore 


r  t 


9a 

9x 


=  2 


J      x,  x2  x3  dt  +  J 


K   x,  x,  dt 
e   2   3 


28 


-  a 


2  2   J  + 

3  X2   dt 


/ 


x3  x2dt 


(30) 


Figure  6  shows  the  analog  implementation  of  equation  (30) . 

As  in  the  identification  of  parameter  "b",  (Section 
III-A),  this  implementation  was  augmented  by  using  the 
properties  of  exponentially-mapped-past  variables  (see 
Section  III-D).   This  was  done  in  order  to  avoid  the  trouble- 
some iterative  reinitialization. 

C.    KALMAN  FILTER 

As  stated  in  Section  II,  the  linear  equations  (15)  and 
(16)  are  used  in  the  Kalman  filter  for  the  estimation  of  the 
state  variable  D.   These  equations  are  repeated  here  for 
convenience . 


X!  =  A(x3  -  Xi) 


(15) 


x2  =  y(x3  "  1  "  bx2)  . 


(16) 


In  the  application  of  the  Kalman  filter  the  variable  x3  is 
treated  as  a  forcing  function. 

Figure  7  shows  the  block  diagram  of  the  estimator  for 
the  system  defined  by  the  difference  equations 

x(k=i)  =  *  x  (k)  +  r  w  (k) 

z(k)  =Hx(k)  +  v  (k) 

where  k  represents  a  point  in  the  discrete  time  domain,  _x(k) 
is  the  vector  of  the  state  variables  at  time  k,  w(k)  is  the 
system  input  which  in  this  case  is  the  sum  of  two  components, 
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the  variable  x3  (a  deterministic  input)  and  Gaussian  noise. 
z.(k)  is  the  vector  of  measured  states  at  time  k,  and  v(k) 
is  the  measurement  noise.   The  matrices  $.  r  and  H  are, 
respectively,  the  state  transition  matrix,  the  distribution 
matrix  and  the  measurement  matrix. 

Recurrence  relations  for  the  optimal  filter,  [14], 
[15] ,  are 

x(k|k)  =  x(k|k-l)+G(k) [z(k) -H  x(k|k-l)] 

x(k+l|k)   =  J  x(k|k), 

where  x(k|k)  is  the  estimated  state  vector  at  time  k  using 
observations  up  to  and  including  time  k,  x(k|k-l)  is  the 
predicted  state  vector  at  time  k,  using  observations  up  to 
time  k-1  and  G(k)  is  the  gain  matrix. 

HI 

The  gain  matrix  G,(k)  is  given  by 

jG(k)  =  ^P(k|k-1)HT[H  P(k|k-1)HT  +  R] "  X 

P(k|k-1)  is  the  covariance  matrix  of  prediction  error  defined 
by 

j?(k|k-l)  =E{[x(k)  -  x(k|k-l)] [x(k)  -  x(k | k- 1) ]T} . 

The  auxiliary  equations  that  determine  the  propagation 
of  P  are 

JP(k|k)  =  (I.  -  £(k)  H)  P(k|k-1) 

^P(k+l|k)   =  £  P(k|k)  j>T  +  Q^ 

where  £(k | k) ,  called  the  covariance  matrix  of  estimation  error, 
is  defined  by 
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P(k|k)  =  E{[x(k)  -  x.(k|k)][x(k)  -  x(k|k)]T}. 


R,  the  covariance  matrix  of  the  measurement  noise,  and  Q, 
the  covariance  matrix  of  the  system  input  noise,  are  de- 
fined as 

R  =  E[v(k)  vT(k)] 

Q  =  f  E[w(k)  wT(k)]£T. 

The  matrix  form  of  the  difference  equations,  that  cor 
respond  to  the  differential  equations  (15)  and  (16) ,  used 
for  the  filter,  is 


x(k+l)  = 


.992    0 


0    .990 


x(k) 


008     0 


02 


J 


*i-\ 


H  = 


1 
u  c 


for  a  sampling  interval  of  .02  sec.  and  a  value  of  .5  for 
the  parameter  "b". 

Since  "bM  is  required  for  the  $  matrix  then  the  iden- 
tification  of  "a"  is  done  after  the  identification  of  b  has 
been  completed. 

The  signal-to-noise  ratio  is  specified  as  the  ratio  of 
two  exponentially-mapped-past  averages,  the  square  of  the 
signal  and  the  square  of  the  noise  signal.   This  average 
function  is  explained  in  the  next  section. 


33 


The  computer  listing  of  the  estimator  algorithm  con- 
taining the  complete  set  of  equations  that  characterize  the 
estimator  is  within  the  subroutine  ESTIM. 

D.    THE  EXPONENTIALLY-MAPPED-PAST  AVERAGE  FUNCTION 

The  open-loop  approach  used  for  the  implementation  of 
the  Kalman  estimator  produced  system  equations  which  do  not 
depend  explicitly  on  the  unknown  parameter  "a".   This  imple 
mentation  of  the  Kalman  estimator  allows  use  of  the  regres- 
sion analysis  technique.   However,  this  procedure  allows 
the  variable  x3  to  reach  the  identifier  without  being  fil- 
tered.  Another  problem  still  present  is  the  iterative  re- 
initialization. 

To  solve  these  problems  an  exponentially-mapped-past 
(EMP)  average  function,  [13],  was  found  very  useful.   The 
EMP  function  also  was  found  to  be  very  suitable  for  acting 
both  as  a  filter  and  eliminating  the  overloading  of  the  in- 
tegrators used  in  the  identification.   This  exponential 
weighting  function  decreases  with  early  values  of  time  and 
its  implementation  allows  for  relatively  fast  control  of 
its  time  constant.   It  gives  the  advantages  of  easy  com- 
putation and  simple  implementation. 

The  EMP  average  function  is  defined  as 

00 

fa(t)  =  a  J   f (t)  ea(T"t)u(t-T)dT 

-00 

where  u(t)  is  a  unit  step  function.   This  convolution  in- 
tegral of  f(t)  and  an  impulse  response  of  a  low-pass  RC 
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filter,  with  leakage  time  constant  RC  =  1/a  is  shown  im- 
plemented in  Figure  8. 

Figures  9  and  10  illustrate  the  changes  made  in  order 
to  include  the  EMP  function  into  the  simulation.   The 
changes  are  the  addition  of  input  potentiometers  with  set- 
tings of  a. 
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f(t) 


R=l/a 
AAW 


Figure  8.   Exponentially-mapped-past  average 
circuit . 
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IV.   EXECUTIVE  ROUTINE 

One  of  the  most  important  tasks  prior  to  obtaining 
data  from  the  analog  simulation  is  the  verification  of  the 
operational  performance  of  the  analog  computer.   This  is 
especially  critical  where  the  analog  computer  used  in  the 
simulation  is  large  and  has  many  other  users,  [20].   This 
task  is  very  difficult  and  no  attempt  was  made  here  to  cover 
the  whole  system;  however,  it  was  found  worthwhile  to  check 
out  all  components  used  in  this  simulation. 

The  executive  routine  MAIN  can  be  divided  into  two 
subsections :   the  first  one  performs  the  checkout  and  the 
second  one  controls  the  actual  identification. 

The  test  performed  by  the  checkout  routines  are: 

1.  TESTO.  -  This  subprogram  tests  the  library  sub- 
routines SETPOT,  SCAN  and  DAC .   To  accomplish  this  TESTO 
sets  potentiometers  and  compares  the  desired  setting  with 
the  actual  setting,  allowing  for  the  inherent  uncertainty 
present  in  the  fourth  digit.   In  the  same  way  the  digital- 
to-analog  converters  (DAC's)  are  tested  for  proper  setting 
and  also  for  any  interaction  among  them.   Any  detected 
error  will  cause  a  message  to  be  printed  and  a  PAUSE  to  be 
executed. 

2.  TEST1.  -  This  subporgram  tests  the  initial  con- 
dition setting  on  the  plant  and  the  analog-digital  conver- 
ters and  digital- analog  converters.   It  takes  1000  samples 
of  each  set  of  converted  values  and  prints  the  mean  errors 
and  the  error  variances  of  both  tests. 
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3.  TEST2.  -  This  subroutine  tests  the  integrators 
used  in  the  simulated  nuclear  reactor  system.   It  prints 
out  the  percentage  error  of  the  integrator  integration 
rates . 

4.  TEST3.  -  This  subprogram  tests  the  initial  condi- 
tions on  the  simulated  plant,  which  allows  verification  of 
the  divider  M003.   It  also  runs  the  plant  to  compare  the 
final  value  with  a  precomputed  numerical  solution.   This 
subroutine  prints  the  percentage  error  for  both  tests. 

5.  TEST4.  -  This  subroutine  performs  a  test  of  both 
identification  systems  using  a  precomputed  numerical  solu- 
tion for  the  steady-state  and  prints  the  percentage  errors 
of  both  systems  compared  with  the  parameter  value  used  in 
the  numerical  solution. 

The  identification  of  the  unknown  parameters  is  per- 
formed by  subroutines  RUN1  and  RUN2. 

1.  RUN1 .  -  Activates  the  identification  system  for 
the  parameter  "b" ,  which  is  performed  first  because  the 
Kalman  filter  requires  the  value  of  the  unknown  parameter 
"b".   The  internal  subroutine  SAMP1,  activated  by  interrupt 
53,  performs  the  sampling  of  x2  and  x3,  computes  the  deriva- 
tive of  x2  by  a  difference  approximation,  sets  the  values 
of  the  reactivity  input  and  the  required  products  in  the 
simulated  nuclear  plant  and  identification  system  respec- 
tively.  RUN1  prints  on  the  teletypewriter  the  value  of  "b" 

2.  RUN2.  -  Activates  the  identification  system  for 
the  parameter  "a".   HANDY,  an  internal  subroutine  activated 
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by  interrupt  53,  samples  X2  and  X3,  calls  ESTIM  (Kalman 
filter)  to  obtain  an  estimate  of  xi ,  performs  the  required 
multiplications  and  sets  these  values  in  the  analog  computer 
for  the  identification  system. 

The  program  was  written  primarily  in  FORTRAN  IV  for 
the  XDS-9300;  some  assembler  language  statements  were  used 
in  order  to  be  able  to  use  the  library  subroutine  DELAY. 

A  sample  of  the  output  taken  under  conditions  of 
failute  of  the  library  hybrid  executive  routine  is  shown 
as  a  part  of  the  Computer  Output. 
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V.   SIMULATION  RESULTS 

Several  tests  were  performed  for  each  identification 
system.   The  first  test  of  each  system  is  shown  in  Fig- 
ure 11  and  Figure  12.   These  figures  show  the  response  of 
the  identification  system  to  a  constant  input  whose  value 
was  numerically  precomputed  and  corresponds  to  a  steady- 
state  value  of  the  plant,  i.e.,  the  digital  multipliers 
shown  in  Figures  9  and  10  were  furnished  with  these  steady- 
state  values.   Since  these  runs  were  performed  without  state 
variable  trajectories  they  do  not  reflect  the  previously 
specified  desired  objective  of  an  identification  system. 
Nonetheless,  for  this  example  the  identification  system 
converges  to  the  actual  parameter  value  if  supplied  with  a 

o 

noiseless  discrete  sample  representing  the  steady-state 
value  of  the  plant  variables.   The  time  scales  are  un- 
labeled since  the  convergence  time  can  be  altered  arbi- 
trarily in  any  noiseless  situation  by  adjusting  the  system 
gain,  without  changing  the  final  result.   Figure  15  shows 
the  type  of  plant  transient  condition  used,  a  truncated 
ramp.   For  the  identification  of  the  unknown  parameter  "b" 
a  low  gain  was  used  and  for  the  parameter  "a"  a  high  gain 
was  used.   These  choices  of  low  and  high  gain  were  made 
arbitrarily  because  for  noiseless  runs  the  procedure  was 
not  sensitive  to  gain  and  either  high  or  low  gains  produced 
the  same  satisfactory  results.   Here  again,  for  the  reason 
stated  previously,  there  is  no  need  for  labeling  the  time 
scales . 


42 


Figures  16  thru  23  show  the  behavior  of  the  system 
under  simulated  real  conditions. 

The  system  response  in  the  identification  of  the  parame 
ter  "b"  with  two  different  levels  of  measurement  noise  is 
shown  in  Figure  16,  where  the  noise  was  added  to  the  vari- 
able that  represents  the  plant  temperature,  x2 . 

Figures  17  and  18  show  the  system  response  in  the  iden- 
tification of  parameter  "a"  with  two  different  levels  of 
measurement  noise.   Unbiased  noise  was  added  to  the  vari- 
able that  represents  the  plant  temperature,  x2 ,  and  to  the 
variable  that  represents  the  neutron  density,  x3. 

Figures  19  thru  21  were  obtained  with  biased  measure- 
ment noise,  and  Figures  22  and  23  show  typical  inputs  to 
the  identification  system. 

Talbe  II   summarizes  the  numerical  results  for  various 
runs  with  different  values  of  the  parameters  "a"  and  "b". 
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TABLE  II 
NUMERICAL  RESULTS 


FAKAM 

a 

b 

0.4 

0.3 

0.4 

0.4 

0.4 

0.5 

0.4 

0.6 

0.4 

0.7 

0.4 

0.8 

0.2 

0.5 

0.3 

0.5 

0.5 

0.5 

0.6 

0.5 

0.7 

0.5 

IDENTIFICATION 


0.40 

0.29 

0.40 

0.39 

0.40 

0.50 

0.41 

0.61 

0.40 

0.71 

0.39 

0.81 

0.21 

0.50 

0.29 

0.50 

0.50 

0.49 

0.59 

0.50 

0.69 

0.50 

The  signal-to-noise  ratios  used  in  all  these  runs  was  20:1 
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TIME+ 


Figure  11. 


Identification  of  parameter  "b"  under  steady 
state  conditions. 


TIME-*- 


Figure  12.   Identification  of  parameter  "a"  under  steady 
state  conditions. 
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conditions . 
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Figure  15.   Reactivity  function  used  in  the  simulation 
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Figure  18.   Identification  of  parameter  "a"  with  measurement 
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Figure  19.   Identification  of  parameter  "b"  under  a  biased 
noise  condition. 
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Figure  22.   Typical  input  signals  for  the  identifier  of 
parameter  "b". 
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Typical  input  signals  for  the  identifier  of 
parameter  "a" . 
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VI.   CONCLUSIONS 

From  the  results  of  the  numerical  example  used,  the 
following  conclusions  are  reached  regarding  the  desired 
goals  for  the  identification  procedure. 

A  typical  response  of  the  identification  system  where 
the  initial  value  of  the  parameters  are  set  to  zero  is  shown 
in  Figures  16  thru  21.   The  results  indicate  that  the  pro- 
cedure does  not  require  a  guessed  value  to  be  effective. 
Other  runs  were  made  with  various  initial  parameter  values 
and  if  the  guessed  value  was  close  to  the  actual  value,  the 
convergence  time  was  greatly  reduced  as  might  be  expected. 
In  any  case,  whether  or  not  a  non  zero  initial  value  was 
used,  the  system  did  not  require  iterative  reinitialization. 

The  simplicity  of  the  identification  procedure  imple- 
mentation, together  with  the  small  number  of  required  analog 
components,  provides  for  a  very  reliable  system,  since  ana- 
log integrators  and  analog  multipliers  are  very  reliable 
components . 

The  average  convergence  time  for  a  signal-to-noise 
ratio  of  10:1  was  found  to  be  approximately  16  seconds,  and 
about  10  seconds  for  noiseless  measurements.   If  the  signal- 
to-noise  ratio  is  small  the  convergence  time  is  large. 

For  the  example  considered  in  this  study  the  procedure 
identified  the  unknown  parameters  within  3%  of  their  actual 
values  with  signal-to-noise  ratios  as  low  as  10:1  in  the 
measured  states. 
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The  identification  system  for  parameter  "b"  was  very 
insensitive  to  the  noise,  except  for  biased  noise  (with  a 
mean  of  101  of  the  stead-state  value  of  the  state  variable) , 
as  shown  in  Figure  20.   The  extreme  biased  condition  was 
taken  as  the  maximum  noise  for  which,  even  when  the  system 
does  not  actually  converge,  the  identifier  output  oscillates 
within  4%  of  the  actual  parameter  value. 

The  identification  of  parameter  "a"  also  showed  very 
good  performance  but  it  was  more  affected  by  measurement 
noise  on  x3  than  on  x2.   This  was  expected  since  x3  is  not 
filtered  by  the  Kalman  filter  and  x3  is  an  important  factor 
in  every  term  of  equation  (30) .   This  problem  can  be  elim- 
inated by  using  filtering  in  the  measuring  circuits  for  x3. 
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APPENDIX  A:   DYNAMIC  EQUATIONS  SIMPLIFICATION 

As  seen  in  the  case  of  an  heterogeneous  reactor  made 
up  of  n  independent  media,  the  state  variables  are  related 
by  the  following  non  linear  differential  equations: 


dn 

at  =  lpe 


Z^j    T.J     jS-   ^n   +    ^A.C.  (1.1) 


dCi        Bi 

— jr  =   7-5-  n    -    X.C.  (1.2) 

dt        A*  1    1  v        J 


■1  ?-4H  1 


•jlCTj    -   Tk).  (1.5) 


The  initial  conditions  for  this  set  of  equations  are 


n(0) 

=    n 
0 

3i 

C.(0) 

=    C       = 
10 

1 

1 

Ji 

0 

T..(0) 

=    0. 

For  the  purpose  of  this  study  the  complexity  of  the 
model  is  reduced  to  a  single  time  constant  approximation, 
i.e.,  i  =  1  and  j  =  1,  so  that 

at  =  (pe  -  aT)  I*  -  r^n  +  AC  f2-« 

||  =   ji  n    -    AC  (2.2) 

■ '&  -  "'fir  -  0   "oT  <2-3> 
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with  the  following  initial  conditions 


Pe(0)  =  0 
n(0)  =  n 


C(0)  =  c 


3n, 


o     l*\ 

T(0)  =  0  . 


A  further  manipulation  will  produce  a  more  convenient 
form,  i.e., 


dn 
dT 


_3 
I* 


fT 


a 


). 


A£* 

+  *£-  c 


B  T  -  l)   n  +  M^  no  C 


or 


1 

dn 

$ 

n 

dt 

I* 

0 

a  T 

3  * 


In  the  same  way 


'       o     o 


(3.1) 


dc 
3t 


y^  n  -  AC 


=  A 


=  A 


(A  n  - c) 


and 


C    dt 
o 


\  o 


•s-^-0  -- 


dt  _  y 
dt    c 


fe-9 


(3.2) 


(3.3) 
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Now   if 


0 

o 

A    U 

1 

A    a 

3  =  i 

then 


N   =   jl      [(Ke-aT-l)    N   +    d] 
D    =    X    (N    -    D) 
T    =    y    (N-l-bT)  , 


and 


N(0)    =    1 
D(0)    =    1 


T(0)    =    0 


Ke(0)    =    0 


A    p 


K 


A   a 


y     ' 


(4.1) 
(4.2) 
(4.3) 
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APPENDIX  B:   REACTOR  STABILITY 

As  it  was  shown  by  Liapunov,  "Problems  general  de  la 
stabilite  du  mouvement , "  the  general  problem  of  stability 
of  motion  relative  to  certain  quantities  can  be  reduced  in 
many  cases  to  a  solution  of  problems  on  stability  of  the 
trivial  solution  of  the  system.   In  theorem  1,  "Direct 
Method,"  he  stated:   "If  the  system  of  differential  equa- 
tions that  defines  the  system  of  differential  equations 
that  defines  the  system  is  such  that  it  is  possible  to  find 
a  positive  definite  function  V,  the  derivative  of  which, 
with  respect  to  time,  calculated  by  virtue  of  the  system 
equations  satisfies  the  inequality  dV/dt  <  0,  then  the  sys- 
tem is  stable." 

The  problem  is  to  investigate  whether  or  not  the 
stationary  state  is  stable  with  respect  to  initial  distur- 
bances . 

The  operation  of  the  reactor  (mathematical  model)  is 
defined  by  equations  (3.1),  (3.2)  and  (3.3),  assuming  that 
the  external  control  is  disconnected  and  neglecting  the 
poisoning  effects  (since  they  represent  a  very  slow  phenom- 
enon) the  reactivity  can  be  considered  as  some  function  of 
temperature,  i.e.: 

6K  =  -  f (t)  . 

First  of  all  it  should  be  noted  that  n  and  C  are  always 
positive  (and  so  are  N  and  D) ,  since  they  cannot  be  negative 
for  physical  reasons.   Now  if  the  initial  conditions  are 
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n  =  n(tQ)  >  0 

C  =  C(t  )  >  0 
o      - 

T  =  T(t  )  , 

then  the  solutions  n(t)  and  C(t)  which  correspond  to  these 
initial  conditions  satisfies  the  inequalities 

n(t)  >  0 

C(t)  >  0 

so  the  stability  will  be  investigated  in  the  domain  D(n  >  0, 
C  >  0) ,  inherent  stability  in  the  presence  of  a  single  group 
of  delayed  neutrons 

^|  =  -f(T)  n  -  jg-  n  +  XC  (5.1) 

§=   *|n  -  AC  (5.2) 

||  =  T(n-l)  -  xbT  (5.3) 

which  obviously  admit  the  special  solution  for  the  stationary 

state 

T  =  Tj   some  final  temperature  (6.1) 

n  =  n   =  1  +  bT  (6.2) 

o 

C  =  C  .  (6.3) 

o 

Now  if  B>0,  £*>0,  X>0,  b>0,  1+b  =0  are  verified  then 
this  solution  for  the  system  is  asymptotically  stable;  for 
physical  reasons  the  first  four  inequalities  of  the  statement 
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are  always  fulfilled  and  the  fifth  is  a  consequence  of  the 
fourth  one . 

Considering  the  Liapunov  function 

T  n  C 

V  = 


f(t)dt  +   J   TQ-nQ)   dn  +  iT^-Co^    dc  (7) 


i  o  o 

in  which  the  constants  T: ,  n   and  C  are  given  by  formulas 

(6.1),  (6.2)  and  (6.3),  the  function  V  becomes  zero  for 

T  =  T,,  n  =  n   and  C  =  C   and  under  the  conditions  given  in 

i '       o  o 

Theorem  (1)  by  M.  V.  Popov  in  [5]  it  is  positive  in  the 
domain  n  >  0 ,  C  >  0 . 

From  equations  (5.1),  (5.2)  and  (5.3)  the  total  deriva- 
tive of  V  with  respect  to  time  is 


dv  (n"nJ  R 

|£  =  f(t)  [r(n-l)-rbT]  +  r  -   °   [-f(T)n-^4  n  +  AC] 

(C  '  Co}     8 
+  r  £-±-    [  J*,  n  -  AC] 


=  -rbTf(T)  -  -|  [  j|n  -  AC][C(n-no)-n(C-Co)] 

rnn      ft  2 

=  -rbTf  (T)  -  _°   ^  n  -  AC  (8) 

and  again  under  the  conditions  of  the  same  theorem,  the 
function  above  is  negative  semi-definite  within  the  range 
n  >  0,  C  >  0  and  become  zero  for 

where  T  -  Ti . 
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Since  dV/dt  is  negative  semi-definite,  it  can  be  con- 
cluded that  the  solution  is  stable  in  the  sense  of  Liapu- 
nov. 
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COMPUTER  OUTPUT 

The  computer  output  shown  here  corresponds  to  different 
runs  where  failures  in  the  XDS-9300  library  hybrid  executive 

package  were  detected. 

ERRr'R  IN  P006  .20 
ERROR  IN  P012  .15 
ERROR  IN  P057        .15 

INTERACTION  BETWEEN  T4?l  AND  T420 
INTERACTION  BETWEEN  T422  AND  T420 
INTERACTION  BETWEEN  T<t23  AND  T420 
INTERACTION  BETWEEN  T424  AND  T420 
INTERACTION  BETWEEN  T425  AND  T420 
INTERACTION  BETWEEN  T426  AND  T420 
INTERACTION  BETWEEN  T^27  AND  T420 
INTERACTION  BETWEEN  T430  AND  T^?0 
INTERACTION  BETWEEN  T<+31  AND  T420 
INTERACTION  BETWEEN  T432  AND  T420 
INTERACTION  BETWEEN  T433  AND  T420 

TEST1 
A/DjD/A  TRUNK  LINES  TEST 
PERCENTAGE  ERRORS: 

.0*398         -.02460 
VARIANCES 

.00000  .00000 


TEST2 
INTEGRATORS  RATE 
PERCENTAGE  ERRORS: 

Til)  1.00012 

T(2)  .00232 


TEST3 
INITIAL  CONDITIONS  PERCENTAGE  ERRORS: 
NEUTRON  CONCENTRATION  .05^93 

DELAYED  NEUTRON  CONCENTRATION       2.00037 
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FINAL  VALUE  PERCENTAGE  ERRORS? 


DELAYED  NEUTRONS  C8NCENTRATI6N 

TEMPERATURE 

NEUTRONS  CONCENTRATION 


.44188 
.28460 
.29163 


TEST4 

IDENTIFICATION  PERCENTAGE  ERRORS! 

PARAMETER  R 

PARAMETER  A 


10.00000 
2.50000 
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U.S.  DEPARTMENT  OF  COMMERCE 
National  Technical  Information  Service 

5285  Port  Royal  Road 
Springfield.  Virginia    22151 


November  29,  1972 


United  States  Naval  Postgraduate  School 
Monterey,  California  93940 


72-5 


Re    AD  738  874 
Gentlemen: 


Solution  of  a  Nuclear  Reactor  Parameter  Identification 
Problem 


We  are  unable  to  provide  the  public  with  satisfactory  copy  of  referenced 
document  for  the  reason(s)  checked  below.   Please  help  us  by  sending 
replacement  copy  together  with  this  letter  by  AIRMAIL  to  address  shown 
below.   If  you  are  unable  to  furnish  replacement,  please  return  this 
form  by  AIRMAIL  indicating  where  the  necessary  replacement  copy  is 
obtainable.   Your  cooperation  will  be  appreciated. 


1.   These  pages  are  missing,  please  furnish:  Page  75. 


I I  2.   These  pages  are  not  reproducible,  please  send  good  pages: 


I I  3.   Document  is  not  reproducible,  please  send  good  c 


opy. 


I I  4.   These  pages  have  been  mutilated,  please  send  replacements: 


□  5.   Other 

Input  Section 

National  Technical  Information  Service 

Department  of  Commerce 

5285  Port  Royal  Road 

Springfield,  Virginia   22151 


Sincerely, 


Chief,  Input  Section 
703-321-8517 


Request  for  Replacement  NTIS-212  (11-72) 


To:  Mr.  Thott 

I  regret  that  our  Library  copy 
is  also  incomplete  in  that  it  lacks 
page  75,  too.   We  have  contacted  the 
advisor  in  a  search  for  the  original 
but  it  appears  that  it  has  either 
been  discarded  or  was  retained  by  the 
author,  who  is  an  officer  in  the 
Peruvian  Navy  and  who  has  now  returned 
to  his  country.    yp 

^u-ia^     A  .  >^c<« 
^Georga  R.  Luckett 
Librarian 
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